
% part 1
fn='Grid_E.out';
fprintf('读取GRDECL')
tic
grdecl  = scan_GRDECL(fn);
toc

% part 2
% 避免双精度带来的误差
grdecl.COORD=round(grdecl.COORD,5);
grdecl.ZCORN=round(grdecl.ZCORN,5);

% part 3
% 确定所有角点坐标
disp('确定角点坐标')
tic
[X,Y,Z,pillarmatrep]=determCORNERPOINTcord(grdecl);
toc

% part 4
% 把标记活动单元和非活动单元的ACTNUM按三维空间顺序重构
if isfield(grdecl, 'ACTNUM')
    actv = reshape(grdecl.ACTNUM ~= 0, grdecl.cartDims);
else
    % No ACTNUM.  Assume all cells active.
    actv = true(grdecl.cartDims);
end
% actv_cp=repelem(actv,2,2,2);

% part 5
% 创建结构体M存储网格维度、活动性、单元和节点信息(此步不考虑单元的活动性)
M.DIM=grdecl.cartDims;
M.ACT=actv;
% 对角点坐标进行唯一化处理
[p,~,ip]=unique([Z(:),Y(:),X(:)],'rows','stable');
% A = [Z(:),Y(:),X(:)];
% p = A(ia,:) 且 A = p(ip,:)。
% ip每一行存储该行节点坐标在唯一化之后的节点序号
p=fliplr(p);
M.POINT.CORD=p;
M.POINT.NUM=size(p,1);

M.ELEME.NUM=prod(M.DIM);

P_INDEX = reshape(ip,size(X));
C_INDEX = reshape(1:M.ELEME.NUM,M.DIM);

[r,c,page] = ind2sub(M.DIM,1:M.ELEME.NUM);
M.ELEME.SUB = [r;c;page].';

fprintf('确定每个单元节点组成，')
tic
M.ELEME.VERTEX = cellfun(@(i,j,k) cpindex2elevert_CNPC(P_INDEX,i,j,k),num2cell(r),num2cell(c),num2cell(page),'UniformOutput',false);
M.ELEME.VERTEX = reshape(M.ELEME.VERTEX,M.DIM);
toc

% part 6
% 删除非活动单元
fprintf('删除无效节点/单元并替换M,P_INDEX,C_INDEX中的节点/单元编号，')
tic
[M,P_INDEX,C_INDEX] = deletinactele_CNPC(M,P_INDEX,C_INDEX);
toc

% part 7
% 定位三棱柱单元
% 检查最大内角是否接近180°
% 如果第一层是三棱柱的话，下面的行都是三棱柱,直接检查柱线顶点是否构成三棱柱
fprintf('定位三棱柱单元,')
tic

pillarxyz=reshape(grdecl.COORD,6,[]);
pillarxyz=pillarxyz';
pillarxy=pillarxyz(:,1:2);
pillarmatrepxy=pillarmatrep(:,:,1);
tpind=zeros(0,3);
zbxy=zeros(0,2);
zbxyind=zeros(0,1);

ordk=[4,1,2,3,4,1];
for i=1:M.DIM(1)
    for j=1:M.DIM(2)
        wz=pillarmatrepxy(2*i-1:2*i,2*j-1:2*j);
        pp=pillarxy(wz([1;2;4;3]),:);
        inag=zeros(4,1);
        for k=2:5
            v1=pp(ordk(k),:)-pp(ordk(k-1),:);
            v2=pp(ordk(k),:)-pp(ordk(k+1),:);
            cs=dot(v1,v2)/norm(v1)/norm(v2);
            inag(k-1)=acosd(cs);
        end
        if ~isempty(find((180-abs(inag))<10, 1))
            tpind=[tpind;i,j];
            lsind=find((180-abs(inag))<10, 1);
            zbxy=[zbxy;pp( lsind ,:)];
            zbxyind=[zbxyind; lsind ];
            % ordk(lsind:lsind+2)
        end
    end
end
[~,ia]=ismember(M.ELEME.SUB(:,1:2),tpind(:,1:2),'rows');
unrgele=find(ia);% unrgele中是三棱柱的编号

toc
%%

% 对每个棱柱单元进行四面体剖分，对非断层单元（四棱柱单元）直接进行剖分
% 对断层单元（三棱柱）单元先分成两个三棱柱，再进行四面体剖分
fprintf('处理第一层单元')
DTlayer=cell(sum(sum(M.ACT(:,:,1))),1);
Faultptlayer=cell(size(tpind,1),1);
tic
for i = 1:M.ELEME.NUM
    if M.ELEME.SUB(i,3)==1
        veti=M.ELEME.VERTEX{i}; % veti为第i个单元的节点编号
        veti=veti'; % 编号为i的单元的节点索引
        pele=M.POINT.CORD(veti,:); % 编号为i的单元的节点坐标
        [ue,iu]=ismember(i,unrgele);
        if ue
            vetils=zeros(7,2);
            vetils(1:4,:)=reshape(veti,[],2);
            vetils(5:7,:)=vetils(1:3,:);

            if mod(iu, size(tpind,1) )==0
                zbi=zbxy(end,:);
            else
                zbi=zbxy(mod(iu, size(tpind,1) ),:);
            end

            [~,ip]=ismember(pele(1:4,1:2),zbi,'rows');
            ipf=find(ip);

            vt1=vetils([ipf,ipf+1,ipf+2],:);
            vt1=reshape(vt1,[],1);
            [~,idvt1]=ismember(vt1,veti);

            vt2=vetils([ipf,ipf+2,ipf+3],:);
            vt2=reshape(vt2,[],1);
            [~,idvt2]=ismember(vt2,veti);

            faultvt=vetils([ipf+3,ipf,ipf+1],:);
            faultvt=reshape(faultvt,[],1);
            [~,ftvt]=ismember(faultvt,veti);
            Faultptlayer{iu}=ftvt;

            pele1=M.POINT.CORD(vt1,:);
            pele2=M.POINT.CORD(vt2,:);

            dt1=delaunayTriangulation(pele1);
            dt2=delaunayTriangulation(pele2);

            DTlayer{i}=[replaceind(dt1.ConnectivityList,idvt1);
                replaceind(dt2.ConnectivityList,idvt2)];

            clear vetils zbi ip ipf vt1 vt2 pele1 pele2 dt1 dt2
        else
            dt=delaunayTriangulation(pele);
            % DT{i}=replaceind(t.ConnectivityList,veti);
            DTlayer{i}=dt.ConnectivityList;
            clear pele dt
        end
    end
end
toc

fprintf('重复第一层的剖分结果')
tic
DT=cell(M.ELEME.NUM,1);
Faultpt=cell(length(unrgele),1);
for i=1:M.DIM(3)
    seleind=sum(sum(M.ACT(:,:,1)))*(i-1)+1;
    teleind=sum(sum(M.ACT(:,:,1)))*i;
    DT(seleind:teleind)=cellfun(@replaceind,DTlayer,M.ELEME.VERTEX(seleind:teleind),'UniformOutput',false);
    sunrgeleind=size(tpind,1)*(i-1)+1;
    tunrgeleind=size(tpind,1)*i;
    Faultpt(sunrgeleind:tunrgeleind)=cellfun(@replaceind,Faultptlayer,M.ELEME.VERTEX(unrgele(sunrgeleind:tunrgeleind)),'UniformOutput',false);
end
toc

Faultpt=unique(cell2mat(Faultpt));
DTmat=cell2mat(DT);
pointdata=zeros(M.POINT.NUM,1);
pointdata(Faultpt)=1;
%%
nodes=M.POINT.CORD;
elements=DTmat-1;
filename='result.vtk';
tic
write_vtk(nodes, elements, filename, pointdata)
toc


